*****************************************************************
* Replication directory for                                   ***
* Prime locations                                             ***
* by Gabriel M. Ahlfeldt, Thilo N.H. Albers, Kristian Behrens ***
* Published in American Economic Review: Insights             ***
*****************************************************************
* 01/2025
* Stata
version 17.0


* This do file generates maps of prime locations by US MSA

* Refined mapping  *************************************************************
	capture mkdir "$temp/GISUS_using"
	
* Loop over CBSAs
	qui u "$data_USMETROS/Raw Numeric Data/METRO LIST/METROS.dta", clear	
	qui tab cbsafp
	local Ncbsafp = r(N) // Total number of cities in list
	levelsof  cbsafp, local(USMETROIDS)	
	local count = 1
	
* Begin loop by metro **********************************************************
	foreach num of local USMETROIDS{	
		di "...Working on CBSA `num', number `count' of `Ncbsafp'..."
		
		* Add merge identifier to PL data set
		u "$temp/grid_PL_output_p${SL}/PL_grids_p${SL}_`num'.dta", clear
			capture drop  mlon
			capture drop mlat
			* Rounded coordinates as identifiers for merging
			gen mlon = round(lon,0.001)
			gen mlat = round(lat,0.001)	
			qui foreach var in cell_id1 cell_id2 cell_id3 id_x id_y { // in case any of these variables are in the data from previous runs
				capture drop `var'
			}
		qui split cell_id, p(_)
		qui destring cell_id2, replace
		qui destring cell_id3, replace
		capture gen double id_x=cell_id2
		capture gen double id_y=cell_id3
		capture drop cell_id1
		capture drop cell_id2
		capture drop cell_id3
		qui save, replace
		
		* Generate data from shapefile
		qui shp2dta using "$temp/NEWGRIDS/GRID_`num'.shp", database("$temp/GISUS_using/db_using_`num'") ///
		coordinates("$temp/GISUS_using/coord_using_`num'") replace // Full grid with employment data
	
	* Prepare the PL grid
		u "$temp/GISUS_using/coord_using_`num'", clear
			qui gen nid =_n // Since merge may affect sorting, we need to save a sorting variable to restore the orginal sorting after merging
			qui sum _X
			local xmin = r(min)
			local xmax = r(max)
			local xrange = r(max)-r(min)
			local width = int(`xrange' / 100000 * 2000) // These adjustments are made to ensure that the maps are plotted in proportionate scale, larger cities have larger maps
			qui sum _Y
			local ymin = r(min)
			local ymax = r(max)
			qui merge m:1 _ID using  "$temp/GISUS_using/db_using_`num'", keepusing( cbsafp id_x id_y) // merge ID variable for mapping mlat mlo
			drop _m 
			qui merge m:1 cbsafp id_x id_y using "$temp/grid_PL_output_p${SL}/PL_grids_p${SL}_`num'.dta"
			sort nid // Restore the orginal sorting
			qui keep if _m == 3
			qui drop _m
			qui drop nid
			qui save "$temp/GISUS_using/PL_coord`num'.dta", replace
			qui egen labX = mean(_X), by(PLID)
			qui egen labY = mean(_Y), by(PLID)
			keep PLID labX labY PLrankByCity  
			qui duplicates drop
			qui save "$temp/GISUS_using/PL_lab`num'.dta", replace		
	
	* Now merge geo ID to data and plot the results
		u "$temp/grid_PL_output_p${SL}/grid_PL_output_p${SL}_`num'.dta", clear
			qui foreach var in cell_id1 cell_id2 cell_id3 id_x id_y { // in case any of these variables are in the data from previous runs
				capture drop `var'
			}
				qui split cell_id, p(_)
				qui destring cell_id2, replace
				qui destring cell_id3, replace
				capture gen double id_x=cell_id2
				capture gen double id_y=cell_id3
				capture drop cell_id1
				capture drop cell_id2
				capture drop cell_id3		
		local title = metro_name[1]
		qui keep cell_id emp developable PLID lat lon  cbsafp id_x id_y
			gen mlon = round(lon,0.001)
			gen mlat = round(lat,0.001)
		qui merge 1:1 cbsafp id_x id_y using  "$temp/GISUS_using/db_using_`num'", keepusing(_ID ) // merge ID variable for mapping
		qui drop _m
		qui merge m:1 PLID using  "$temp/GISUS_using/PL_lab`num'.dta" 
		qui drop _m
		qui replace emp = 0 if emp == .
		qui replace emp = -1 if developable == 0
		qui replace labX= labX+0.015
		qui replace labY= labY+0.015
		qui global PCT -1.5 -0.5 0.5 1.5 50 100 250 1000 1000000 // defines the cutoffs in the legend
		qui spmap emp using "$temp/GISUS_using/coord_using_`num'", id(_ID) fcolor(gs14 gs12 gs10 gs8 gs6 gs5 gs4 gs2 ) ocolor(white ..) osize(0.00000 ..)   ///
				clmethod(custom)  clbreaks($PCT)  legend(off) ///
				 label(xcoord(labX) ycoord(labY) label(PLrankByCity)   color(red) size(4) )  ///
				 legend(size(3) region(fcolor(gs16))) 	///
				polygon(data("$temp/GISUS_using/PL_coord`num'") fcolor(maroon%75)  ocolor(maroon)  osize(0.00000 ..) ) 
* Export 2D maps of prime locations for all MSAs, including those in Figure 1
					capture mkdir "$figures_App/US-CBSAs"
					capture mkdir "$figures_App/US-CBSAs/"
					capture mkdir "$figures_App/US-CBSAs/GISMAPS"
					qui graph export "$figures_App/US-CBSAs/GISMAPS/MAPS_PL_TotalEmp_BaseP_`num'.png", replace height(1000) width(1000)	//  width(`width')  width(1000)		
		capture erase "$temp/GISUS_using/db_using_`num'.dta"
		capture erase "$temp/GISUS_using/PL_coord`num'.dta"
		capture erase "$temp/GISUS_using/PL_lab`num'.dta"
	* Update counter
		local count = `count'+1
	}
* LOOP ENDS ********************************************************************

// Gen folder for saving selected maps for Figure 1
	capture mkdir "$figures/US-CBSAs/Figure_1" 

	// Copy additional file to Figure_1 folder
	capture mkdir "$figures/US-CBSAs/Figure_1" 
	copy "$figures_App/US-CBSAs/GISMAPS/MAPS_PL_TotalEmp_BaseP_16980.png" /// 	
	"$figures/US-CBSAs/Figure_1/MAPS_PL_TotalEmp_BaseP_16980.png", replace	
	copy "$figures_App/US-CBSAs/GISMAPS/MAPS_PL_TotalEmp_BaseP_35620.png" /// 		
	"$figures/US-CBSAs/Figure_1/MAPS_PL_TotalEmp_BaseP_35620.png", replace	
	copy "$figures_App/US-CBSAs/GISMAPS/MAPS_PL_TotalEmp_BaseP_31080.png" ///   
	"$figures/US-CBSAs/Figure_1/MAPS_PL_TotalEmp_BaseP_31080.png", replace	

* Script ends